module smooth_indicators_mod

  implicit none

  interface
    pure real(8) function smooth_indicator_interface(a)
      real(8), intent(in) :: a(:)
    end function
  end interface

contains

  pure real(8) function smooth_indicator_3(a) result(res)

    real(8), intent(in) :: a(:) ! 3

    res = a(2) * a(2) + 13 * a(3) * a(3) / 3

  end function smooth_indicator_3

  pure real(8) function smooth_indicator_2x2(a) result(res)

    real(8), intent(in) :: a(:) ! 4

    res = a(2)**2 + a(3)**2 + 7 * a(4)**2 / 6

  end function smooth_indicator_2x2

  pure real(8) function smooth_indicator_3x3(a) result(res)

    real(8), intent(in) :: a(:) ! 9

    res = (  720 * a(2) * a(2) + 3120 * a(3) * a(3) + 720  * a(4) * a(4) + 840   * a(5) * a(5) &
          +  120 * a(4) * a(6) + 3389 * a(6) * a(6) + 3120 * a(7) * a(7) + 120   * a(2) * a(8) &
          + 3389 * a(8) * a(8) + 520  * a(3) * a(9) + 520  * a(7) * a(9) + 13598 * a(9) * a(9) ) / 720

  end function smooth_indicator_3x3

  pure real(8) function smooth_indicator_4x4(a) result(res)

    real(8), intent(in) :: a(:) ! 16

    res = a(2)**2 + (13*a(3)**2)/3 + (3129*a(4)**2)/80 + a(5)**2 + (7*a(6)**2)/6 + a(5)*a(7)/6 + (3389*a(7)**2)/720        &
          + (17./30)*a(6)*a(8) + (47459*a(8)**2)/1120 + (13*a(9)**2)/3 + a(4)*a(10)/24 + (3389*a(10)**2)/720               &
          + (13./18)*a(3)*a(11) + (13./18)*a(9)*a(11) + (6799*a(11)**2)/360 + (1043./160)*a(4)*a(12) + (73./32)*a(10)*a(12)&
          + (22846129*a(12)**2)/134400 + a(2)*(12*a(4) + 4*a(10)/24 + a(12)) + 0.5*a(5)*a(13) + a(7)*a(13)/24              &
          + (3129*a(13)**2)/80 + (17./30)*a(6)*a(14) + (11./80)*a(8)*a(14) + (47459*a(14)**2)/1120 + a(5)*a(15)/24         &
          + (73./32)*a(7)*a(15) + (1043./160)*a(13)*a(15) + (22846129*a(15)**2)/134400 + (11./80)*a(6)*a(16)               &
          + (114997*a(8)*a(16))/5600 + (114997*a(14)*a(16))/5600 + (19583517*a(16)**2)/12800
    
  end function smooth_indicator_4x4

end module smooth_indicators_mod
